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ABSTRACT 


Data  from  simulations  of  channel  flow  and  decay  of  homogeneous  turbulence  in¬ 
dicate  anomalously  strong  correlation  of  velocity  and  vorticity  directions  ('local  Beltram Na¬ 
tion ’)  in  band-filtered  velocity  fields  when  the  band  consists  of  a  thin  cigar  in  mode  space 
(whose  physical  space  representation  is  as  an  array  of  'pancake  eddies’).  Spherical  shells  or 
other  broad  bands  in  mode  space  do  not  seem  to  exhibit  the  effect. 


In  fully  developed  three-dimensional  turbulence,  ’local’  interactions  (that  is,  in¬ 
teractions  between  comparable  scales)  have  been  traditionally  thought  to  be  dynamically 
much  more  significant  than  nonlocal  interactions  (for  example,  in  the  cascade  of  turbulent 
energy).1 

However,  analysis  of  experimental  data  on  energy  transfer  shows2  that  those  in¬ 
teractions  which  contribute  most  strongly  involve  wavenumber  triads  which,  while  local, 
have  aspect  ratios  (ratio  of  largest  to  smallest  wavenumber)  in  the  range  5  to  10.  The 
ALHDIA3  yields  results  typical  of  second-order  turbulence.  It  predicts  that  about  20%  of  the 
total  transfer  come  from  triads  interactions  with  aspect  ratios  less  than  2,  while  about  80% 
come  from  triads  with  aspect  ratios  less  than  eight.  However,  these  figures  give  less  than  the 
whole  story.  If  local  interactions  are  consistently  removed  from  the  dynamics,  the  loss  of  the 
energy  transfer  from  these  interactions  is  counterbalanced  with  an  enhanced  transfer  by  the 
remaining  interactions,  at  least  within  the  framework  of  the  closures5.  The  result  is  that  re¬ 
moving  triads  with  aspect  ratios  eight  and  less  leaves  the  total  energy  transfer  nearly  un¬ 
changed,  instead  of  reducing  it  by  80%. 

To  analyze  the  effect  of  local  versus  nonlocal  interactions,  we  decompose  the 
velocity  field  v(r)  and  the  vorticity  field  o>  =  V  x  v  as  v(r)  =  £v>(r);  oo(r)  =  • 

Here  the  band-filtered  velocity  field  vJ(r)  is  defined  as  v'(r)  =  £u(k)exp(ik-r)  where  the 
sum  extends  over  wavcvcctors  k  within  suitable  distinct  waveveetor  bands  B;.  More  general 
ly,  the  band  filtered  field  is  defined  as  vJ(r)  =  V  u(n)0n(r)  where  <>n(r)  is  a  complete  set  <>t 


orthogonal  functions.  In  the  channel  flow  simulations  reported  below,  the  velocity  field  is 
Fourier  transformed  in  the  x-  and  y-directions  while  the  Chebyshev  polynomials  Tp(z)  are 
used  in  the  z-direction  perpendicular  to  the  walls.  In  the  simulations  of  decaying  turbulence 
reported  below,  three-dimensional  Fourier  series  are  used. 

In  terms  of  the  band-filtered  fields,  the  Navier-Stokes  equations  for  incompressi¬ 
ble  flow  are 

4f«ZNy+V„V*V  (,) 

i.j 

where  Ny  =  v‘  x  ©i  -  V  n y  and  V  tty  subtracts  out  the  divergent  part  of  v1  x  to’.  A  meas¬ 
ure  of  the  strength  of  local  interactions  is  given  by  the  intra-band  contribution  Nj  {  to  the 
nonlinear  term  in  (1). 

In  two-dimensional  flow,  velocity  and  vorticity  are  perpendicular,  so  the  intra- 
band  term  Ny  can  only  be  small  through  pressure  balancing.  This  is  indeed  observed  in 
high-Reynolds  number  numerical  simulations  of  decaying  two-dimensional  flow,  which 
display  organization  of  the  flow  into  quasi-one-dimensional  vortex-gradient  sheets8  or  nearly 
circular  vortex  patches9. 

In  three-dimensional  flow,  depletion  of  local  nonlinear  interactions  can  occur  ei¬ 
ther  by  pressure  balance  within  N( ,  or,  more  directly,  by  local  ’Beltramization '  in  which  the 
band-filtered  fields  vl  and  a)1  are  nearly  parallel  in  physical  space.  While  local  pressure  bal 


ance  does  not  allow  for  complicated  instantaneous  three-dimensional  flow  topologies, 111  local 
Beltramized  flows  are  not  so  restricted.  In  this  work,  we  have  investigated  the  possibility  of 
local  Beltramization  using  results  of  a  direct  numerical  of  turbulent  channel  flow  at  a  Rey¬ 
nolds  number  Re*  =  184  based  on  the  half-channel-width  and  the  friction  velocity.  Thus  the 
channel  width  in  dimensionless  wall  coordinates  is  H  =  368.  The  computer  code  is  described 
in  Ref.  11.  Two  channels  of  different  linear  dimension  have  been  considered.  The  number 
of  grid  points  used  in  these  simulations  was  323  and  322x  64,  respectively.  The  same  effect 
has  also  been  studied  in  (1283)  simulations  of  the  decay  of  homogeneous  turbulence. 

The  analysis  used  here  follows  an  earlier  analysis  of  the  full  (non-bandpass- 
filtered)  flow  field  v(r).12  Specifically,  we  evaluate  the  distribution  of  the  cosine  of  the  angle 
0  between  the  velocity  v'(r)  and  vorticity  o)’(r).  We  subdivide  the  interval  - 1  <;  cos0  <  1 
into  two  hundred  equally  spaced  intervals  and  define  the  probability  function  P(cos0)  as  the 
number  of  grid  points  having  the  angle  0  between  the  bandpass-filtered  velocity  and  vortici¬ 
ty.  The  variable  cos  0  is  the  relevant  one  in  three-dimensions  because  of  the  form  of  the 
volume  element.  Note  that  £P(cos0)  =  N  where  N  is  the  total  number  of  grid  points.  Each 
half  of  the  channel  (0  <  z  <  184)  is  subdivided  into  three  intervals:  0  <  z  <  15: 
15  ^  z  <  40;  40  <  z  <,  184.  In  order  to  avoid  the  effects  of  the  walls,  the  velocity  field  has 
been  analyzed  in  the  outer  part  of  the  channel  40  <  z  <  328.  The  Fourier-Chebyshev 
band-filtering  has  been  performed  using  the  data  on  the  entire  flow  field.  Then  the  filtered 
field  is  transformed  back  into  physical  space  and  the  resulting  ticld  is  analysed  locally  in  space 


for  40  <  z  <  328.  We  recognize  that  there  may  be  ambiguity  in  space  localization  due  to  an 
"uncertainty  principle". 


The  results  plotted  in  Fig.  1  for  the  band  B;  consisting  of  kx  =  ±  4;  ky  =  ±  2; 
p=  15,16,17,18  show  that  v'(r)  and  co'(r),  tend  to  align;  that  is  the  probability  function 
P(cos9)  is  sharply  peaked  at  cos0  =  - 1  at  t  =  50.  To  characterize  the  effect  we  define  the 
Beltrami  ratio  as  B=  max[P(l),  P(- 1)]/Pav(0),  where  Pav(0)  is  the  value  of  P(cos0)  aver¬ 


aged  over  the  interval  -  — -  <  cos0  <  — In  Fig.  2,  we  plot  the  time  evolution  of  B(t)  for 


40  <  t  <  77  with  samples  taken  at  the  time  intervals  At  =  1.  We  find  that  B  >  12  approxi¬ 
mately  30%  of  the  time.  Continuing  to  add  modes  to  the  band  leads  to  disappearance  of  the 
peak  of  P(cos0)  when  the  band  is  too  wide.  A  similar  effect  has  been  observed  in  numerical 
simulation  of  decaying  homogeneous  turbulence  (see  Fig.  3).  It  is  typical  for  a  variety  of 
realizations  of  band-filtered  velocity  fields  in  both  channel  and  decaying  turbulence  simula¬ 
tions  that  as  soon  as  a  few  modes  are  present,  there  is  an  organization  of  the  flow  field  in 
which  velocity  and  vorticity  are  almost  aligned. 


It  should  be  mentioned  that  the  band-filtered  velocity  field  v*  does  not  satisfy  the 
incompressibility  condition  V  •  =  0  because  of  the  properties  of  the  Chebyshev  polynomi¬ 

als.  To  assess  the  role  of  incompressibility  we  have  subtracted  the  nonsolenoidal  part  trout 
the  filtered  velocity  field  and  calculated  the  probability  function  P(cos0)  of  the  residual  lick! 


satisfying  the  V  v  =  0  constraint.  The  results  are  presented  in  Fig.  1(b).  2(b).  It  is  <>1hi- 


ous  that  imposition  of  the  incompressibility  condition  docs  not  strongly  change  the  results 


presented  in  Fig.  1(a),  2(a)  since  the  filtered  velocity  field  v(kx,ky,p)  is  only  weakh 
compressible  in  the  part  of  the  channel  we  consider  here  (z+  >  40).  It  is  only  in  the  wall  re¬ 
gion  z+  <  15  that  the  effect  is  not  observed.  It  is  clear  from  Fig.  2(b)  that  imposing  the  in¬ 
compressibility  constraint  leads  to  an  increase  of  the  value  of  B  by  some  10-40%.  An  in¬ 
teresting  property  of  the  probability  density  P(cos0)  is  its  asymmetry  P(l)  *  P(—  1 )  when 
the  Beltrami  ratio  B  is  large.  However,  the  mean  distribution  P(cos0)  obtained  from  35 
realizations  in  the  interval  40  <,  t  <  75  is  very  symmetric. 

This  ’local  Beltramization’  effect  is  observed  to  occur  for  general  wavevector 
packets  W  with  large  aspect  ratio  k  >>  Ak  >>  1  where  k  is  a  typical  wavenumber  in  the 
band  and  Ak  is  the  ’diameter’  of  the  packet  [Ak  =  max(  |k  -  p|  with  k,pe  W)].  These 
packets  produce  ’pancake’  eddies  in  which  both  the  velocity  and  vorticity  are  quasi-two  di¬ 
mensional  fields.  This  may  be  most  easily  seen  by  considering  a  prototype  of  such  a  packet, 
namely,  one  for  which  kz<<  kx,ky  for  all  k  in  the  packet  and  kz  restricted  to  a  narrow 
range.  Thus  the  incompressibility  condition  gives  w  =  -(kxu  +  kyv)/kz  <<  |u  |  +  |v  j,  while 
the  vorticity  field  co  =  ikxv  satisfies  coz  <<  |cox|  +  |coy|.  For  such  pancake  eddies,  we  would 
expect,  in  the  absence  of  other  correlations,  that  the  angle  between  v  and  co  is  uniformly  dis- 


5  tributed  in  the  plane  of  the  pancake,  so  that  P(c os  0)  «  3 /  j s i n  0 !  with  strong  peaks  at 


cos  0  =  ±  1 . 


To  test  whether  the  observed  effect  (sec  Fig.  1,  2)  is  of  dynamic  or  kinematic  na 
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ture  we  chose  a  number  of  distinct  band-filtered  fields  corresponding  to  a  few  different  in¬ 
stants  of  time  t  (see  Fig.  2(a)).  From  these  fields  we  generated  1000  different  realizations  by 
uniformly  distributing  the  complex  numbers  v(kx,ky,p)  for  each  selected  mode  in  the  band. 
Imposing  the  incompressibility  condition  on  each  of  the  band-filtered  random  fields  we  have 
generated  another  1000  realizations.  An  additional  500  random  realizations  were  generated 
by  uniformly  distributing  the  random  complex  numbers  to  each  of  the  (32x32x64)  modes 
of  the  total  field  v(kx,ky,p).  The  resulting  random  velocity  fields  were  made  incompressible 
and  then  band  filtered.  The  weakly  compressible  band-filtered  fields  were  again  made  in¬ 
compressible  by  subtraction  of  the  nonsolenoidal  part.  Then  the  probability  density  P(cos0) 
was  calculated  for  each  member  of  the  ensemble  consisting  of  2500  realizations.  We  then 
calculated  the  cumulative  probability  p(B)  that  gives  the  probability  that  one  Beltrami  ratio  is 
larger  than  B.  The  function  p(B)  derived  from  these  2500  fields  was  compared  with  the 
function  p(B)  corresponding  to  35  realizations  obtained  from  the  time  evolution  of  the 
band-filtered  numerical  solution  of  the  Navier-Stokes  equation  at  40  £  t  <  75  taken  at 
At  -  1.  The  results  are  presented  in  Fig.  4.  These  runs  suggest  that  the  observed  effect 
does  not  have  a  simple  kinematic  explanation.  However  we  must  note  that  the  Chebyshev- 
Fourier  modes  do  not  diagonalize  the  two-point  covariance  function  of  the  observed  velocity 
field  and  for  this  reason  the  randomized  comparison  fields  do  not  exhibit  the  two-point  co- 
variance  of  the  observed  fields.  This  property  will  be  investigated  later13. 


The  generality  of  the  above  results  on  the  bandpass-filtered  geometrical  order  ol 
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velocity  and  vorticity  is  yet  to  be  tested  on  other  turbulent  flows.  It  is  difficult  now  to  give  a 
systematic  explanation  for  local  Beltramization.  This  may  actually  be  as,  or  more,  difficult 
than  explaining  the  remarkable  buildup  of  correlations  between  velocity  and  magnetic  fields 
in  strong  MHD  turbulence.14  We  stress  that  it  is  unlikely  that  any  of  the  traditional  two- 
point  closures  can  provide  insight  into  the  process  of  local  Beltramization.15 

We  now  speculate  on  several  instability  mechanisms  that  may  be  important  in  the 
dynamics  that  lead  to  local  Beltramization.19  First,  secondary  instability  is  a  short 
wavelength  instability  of  quasi-two-dimensional  flows  that  tends  to  produce  locally  aligned 
vorticity  and  velocity.  On  the  other  hand,  anisotropic  pancake  eddies  are  subject  to  the 
long- wavelength  hyperscale  instability.  These  eddies  are  also  unstable  to  the  AKA  instabili¬ 
ty18  which  also  leads  to  generation  of  Beltrami  flows  at  large  scales.  Perhaps  a  combination 
of  these  instabilities,  acting  cyclically,  can  explain  the  observed  effects. 

The  importance  of  steady  solutions  of  the  Euler  equation  in  the  context  of  tur¬ 
bulence  has  been  discussed  by  H.K.  Moffat20,  he  argues  that  turbulent  flows  can  be 
represented  in  terms  of  coherent  structures  with  large  helicity  (Beltrami  flows).20  Our 
findings  support  Moffatt’s  conjecture  in  the  sense  that  the  flow  may  be  a  superposition  of 
Beltrami-like  structures  defined  at  many  different  scales. 

Whatever  the  cause  of  local  Beltramization.  it  can  have  far  reaching  conse¬ 
quences.  Indeed,  we  may  conjecture  that  fulls  developed  turbulence  rearranges  itself  into  a 
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hierarchy  of  coherent  near-Beltrami  flows  with  minimal  self-interaction.  In  this  case,  tur¬ 


bulence  should  be  amenable  to  multiple-scale  perturbation  theory21  or  renormalization-group 
methods.  Furthermore,  virtually  all  turbulence  models,  starting  with  those  of  Prandtl,  are 
based  on  scale  separation  to  justify  eddy  viscosity  concepts.  The  present  ideas  seem  to  pro¬ 
vide  some  justification  for  these  approaches. 
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ficure  captions 


Fig.  1 


Fig.  2 


Fig.  3 


Fig.  4 


Probability  distribution  P(cos0)  in  the  channel  flow  in  the  region  40  <  z+  <  328 
for: 

a.  the  band-filtered  velocity  field  (kx  =  ±  2;ky  =  ±  2;P  =  15,16,17,18) 

b.  the  same  band-filtered  field  with  the  incompressibility  condition  V  •  \i  =  0  im¬ 
posed. 

Time  evolution  of  the  Beltrami  ratio  for  the  band-filtered  field: 

Solid  line:  velocity  field  vj(r),  (kx  =  ±  2;  ky  =  ±  2;  P  =  15,16,17,18) 

Dotted  line:  the  same  field  with  the  incompressibility  condition  imposed. 

Probability  function  P(cos8)  for  the  band-filtered  field  v^ 

(kx  =  ±  2;ky  =  ±  2;k2  =  ±  22)  in  a  simulation  of  decaying,  homogeneous  tur¬ 
bulence. 

Cumulative  probabilities  p(B): 

Curve  I  (•  ):  corresponds  to  31  realizations  of  v(r)  (kk  =  ±  2;ky  =  ±  2;p  =  15-18) 
taken  from  the  time  evolution  of  the  flow. 

Curve  II  (• ):  p(B)  generated  by  introducing  random  phases  into  the  fields  vj(r). 
Curves  (x):  same  as  (•  )  but  with  the  incompressibility  imposed  on  all  fields  vJ(r). 
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It  is  shown  that  the  renormalization  group  theory  of  turbulence  leads  to  the  relation 

Ba  =  CKP,  between  the  turbulent  Prandtl  number  P,,  the  Kolmogorov  constant  CK,  and  the 

Batchelor  constant  Ba. 


The  purpose  of  this  Letter  is  to  show  that  the  renormal¬ 
ization  group  theory1,2  of  turbulence  leads  to  a  simple  rela¬ 
tion  between  the  Batchelor  (Ba)  and  Kolmogorov  (CK) 
constants.  Here  Ba  and  CK  are  defined  by  the  inertial  range 
kinetic  energy  and  passive  scalar  spectra  [£(&)  and  Er(k), 
respectively  ] : 

E  =  CK€2lik  ~sn  (1) 

and 

Et  —  Ba(.V/e,/3)£  -5/3.  (2) 

The  parameters  A'  and  ?  are  the  scalar  and  kinetic  energy 
dissipation  rates  defined  as 

—dr 


2 


Sl+IslY 

dx, J 


and  k0  and  v0  are  molecular  diffusivity  and  viscosity,  respec¬ 
tively. 

In  the  first  step  of  the  renormalization  group  procedure, 
we  eliminate  modes  with  wavenumbers  larger  than  p  from 
the  equation  of  motion  for  the  Founer  components  of  the 
velocity  and  scalar  fields  tz ,  ( kua> )  and  TTk.a),  respectively. 
In  this  calculation,  it  is  assumed  that  k  4p.  In  the  limit  k  4p. 
the  sole  effect  of  the  small-scale  elimination  procedure  is  the 
generation  of  a  tubulent  viscostty  v( p )  and  diffusivity  Kip ) . 
In  the  limit  k-t,pK.kd,  where  kd  is  the  Kolmogorov  dissipa¬ 
tion  wavenumber,  the  turbulent  transport  coefficients  are 
proportional  and 

vlp)  =  P,K{p)  (5) 

with  P,  —  0.7179. 

We  define  the  velocity  and  scalar  fields  vf  and  Tp  fol¬ 
lowing  the  elimination  of  wavenumbers  larger  than  p.  It  can 
be  shown1 2  that  the  two  scalars 


,  vipifthf.  dif 

e(p)  =  — —  —  *  — 
2  \Sx  dx 


PHys  =  uifls  30  ft),  .anoar,  ' 98T 


S(p)  =  Kip) 


£Tf\2 

.  dx,  / 


are  constants  independent  of  p.  This  has  some  important 
implications.  If  the  spectra  are  given  by  relations  ( 1 )  and 
(2)  then 

?  —  2 v(p)  P k  zE(k)dk  =  — v(p)CKe2lipin  (8) 

Jo  2 

and 

.V=2*(/>)f  k:ET{k)dk  =  ^-K(p)^f—Ba.piJ3.  (9) 

Jo  2 

Dividing  (8)  by  (9)  and  using  relation  (5)  we  obtain 

Ba  —  CKP,.  (10) 

Substituting  the  values  CK  =  1.6 17  and/5,  =  0.7179  derived 
in  Refs.  1  and  2  we  obtain  Ba  =  1.161  in  good  agreement 
with  experimental  data. 

This  result  is  obtained  using  the  e-expansion  in  the  vi¬ 
cinity  of  the  fixed  point.  Thus,  relation  (9),  is,  strictly  speak¬ 
ing,  valid  in  the  limit  of  Reynolds  number  R  —  w .  We  believe 
that  formula  (9)  is  accurate  when  molecular  diffusion  is 
negligible  in  comparison  with  turbulent  diffusion.  It  should 
be  mentioned  that  experimental  data  on  Ba  are  much  more 
scattered  than  data  on  the  Kolmogorov  constant  CK .  This 
can  be  easily  explained  since,  in  many  cases,  the  contribution 
from  even  weak  natural  convection  can  strongly  influence 
the  heat  transfer  but  not  a  momentum  transfer  in  turbulent 
flow  when  R  is  not  sufficiently  large. 
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ABSTRACT 


An  efficient  method  for  computing  a  given  number  of  leading 
eigenvalues  (i.e.  having  largest  real  parts)  and  the  corresponding  eigen¬ 
vectors  of  a  large  asymmetric  matrix  M  is  presented.  The  method  con¬ 
sists  of  three  main  steps.  The  first  is  a  filtering  process  in  which  the 
equation  x  -=  Mx  is  solved  for  an  arbitrary  initial  condition  x(0)  yielding: 
x(t)  =  eMtx(0).  The  second  step  is  the  construction  of  (n  +  1)  linearly 
independent  vectors  vm  =  Mmx,  0  <,  m  £  n  or  vm  =  emMtx  (x  being  a 
"short"  time  interval).  By  construction,  the  vectors  vm  are  combinations 
of  only  a  small  number  of  leading  eigenvectors  of  M.  The  third  step 
consists  of  an  analysis  of  the  vectors  {vm}  that  yields  the  eigenvalues  and 
eigenvectors. 

The  proposed  method  has  been  successfully  tested  on  several 
systems.  Here  we  present  results  pertaining  to  the  Orr-Sommerfeld 
equation.  The  method  should  be  useful  for  many  computations  in  which 
present  methods  are  too  slow  or  necessitate  excessive  memory.  In  par¬ 
ticular,  we  believe  it  is  well  suited  for  hydrodynamic  and  mechanical  sta¬ 
bility  investigations. 
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I.  INTRODUCTION 


Numerous  problems  in  science  and  technology  involve  the  computation  of 
leading  eigenvalues  (and  sometimes  of  the  corresponding  eigenvectors)  of  large  asym¬ 
metric  matrices.  In  some  cases,  such  as  that  of  the  calculation  of  eigenvalues  of 
transfer  matrices  in  lattice  problems  (Stanley,  1971)  or  in  hydrodynamic  stability  ana¬ 
lyses  (Drazin  &  Reid,  1981),  the  matrix  of  interest  is  infinite,  in  principle. 

While  rather  efficient  methods  exist  for  the  diagonalization  of  symmetric 
matrices  (see,  e.g.  Parlett,  1980  or  Lewis,  1977)  —  no  satisfactory  algorithm  for  the 
asymmetric  case  is  known  to  us.  The  main  obstacle  is,  of  course,  the  nonorthogonali¬ 
ty  of  the  eigenvectors.  The  Amoldi  (1951)  method  and  its  variants  enable  the  calcu¬ 
lation  of  some  eigenvalues,  but  not  necessarily  those  with  the  largest  real  parts  (which 
are  important  in  stability  investigations;  see  however  Jenning  &  Stewart,  1975  and 
Saad,  1980  and  references  therein).  At  best  this  method  produces  eigenvalues  of  larg¬ 
est  moduli.  Another  disadvantage  of  the  Amoldi  method  is  the  necessity  to  increase 
the  dimension  of  the  Krylov  (Wilkinson,  1965)  subspace  considered  in  order  to  im¬ 
prove  accuracy.  This  may  create  memory  problems. 

Another  method  (Bennetin  et  a l,  1980;  Shimada  &  Nagashima,  1979), 
which  was  designed  to  compute  Liapunov  exponents  for  dynamical  systems,  is  capable 
of  producing  the  real  parts  of  the  leading  eigenvalues.  Its  advantage  is  in  the  fact  that 
its  implementation  does  not  require  a  large  memory.  However  the  method  is  slowly 
converging  and  produces  neither  the  imaginary  parts  of  the  eigenvalues  nor  the  eigen¬ 
vectors.  (However,  see  Goldhirsch  et  a l,  1987,  for  an  accelerated  convergence 
method  and  computation  of  eigenvectors.  The  latter  can  be  computed  efficiently  tor 
relatively  small  systems.) 


It  thus  seems  that  many  important  problems  involving  'urge  asymmetric 
matrices  are  very  difficult  to  analyze  using  available  methods.  The  algorithm  proposed 
in  this  paper  is  designed  for  such  problems.  It  is  fast  and  simple  and  can  therefore  be 


easily  implemented. 

The  proposed  method  has  three  essential  steps.  The  first  step  is  designed 
to  filter  out  nonleading  eigenvectors.  Let  x  be  an  "initial"  vector  and  M  the  matrix 
whose  leading  eigenvalues  are  sought.  We  first  solve  the  equation: 

41=Mx  (1.1) 

for  0  <,  t  <,  t0.  The  resulting  vector  x(t0)  =  eMt°x  is  obviously  a  combination  of 
eigenvectors  corresponding  to  leading  eigenvalues  (henceforth  called  leading  eigenvec¬ 
tors).  The  nonleading  eigenvectors  are  being  damped  by  exponential  factors.  More¬ 
over,  if  it  so  happens  that  the  chosen  initial  vector  x  is  independent  of  a  leading 
eigenvector,  then  this  eigenvector  will  be  introduced  into  x(t)  by  roundoff  errors  in 
the  process  of  numerical  integration  of  eqn.  (1.1).  The  next  step  involves  the  creation 
of  n  linearly  independent  vectors.  This  can  be  done  either  by  computing  (Mmx(t0); 
0  ^  m  <  n-l}orby  computing  (emMTx(t0);  0  £  m  £  n-  1).  It  is  obvious  (and 
shown  below)  that  the  n  vectors  created  this  way  are  independent  for  nondegenerate 
spectra.  A  method  for  degenerate  spectra  will  be  described  below  as  well.  Once  w-e 
obtain  n  independent  vectors  which  are  essentially  linear  combinations  of  leading 
eigenvectors  only,  it  is  a  matter  of  straightforward  algebra  to  produce  the  leading 
eigenvalues  and  eigenvectors. 

The  structure  of  the  paper  is  as  follows.  Section  II  presents  a  description 
of  the  me;hod  we  propose.  Section  III  offers  error  estimates  and  analyses  ol  the  alge- 
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braic  structure  of  the  vectors  produced  in  the  second  step  of  the  method.  Section  IV 
presents  two  algorithms  based  on  the  proposed  method.  Section  V  presents  results 
obtained  from  an  implementation  of  the  method  in  the  case  of  the  Orr-Sommerfeld 
equation.  Section  VI  offers  a  brief  summary.  Those  readers  not  interested  in  the  de¬ 
tails  of  the  mathematical  analysis  of  the  method  may  skip  Sec.  III. 


V.V.V.V.V.  .\ 


II.  DESCRIPTION  OF  THE  METHOD 


Consider  a  matrix  M,  which  is,  in  general,  asymmetric  and  large  (say  of 
rank  R).  Assume  that  the  eigenvalues  of  M:  X1(  X2,  X3  •  •  •  are  arranged  so  that 
Re^  >  ReX2  >  ReX3  •  •  •  ,  i.e.  it  is  assumed  that  the  real  part  of  the  spectrum  of  M 
is  nondegenerate.  Let  e3,  e2,  e3  •  •  •  be  the  right  eigenvectors  of  M  corresponding  to 
the  eigenvalues  Xj,  X2,  X3  •  •  •  ,  respectively. 


Any  vector  x  can  be  expanded  in  terms  of  the  right  eigenvectors  of  M: 


x  =  X  aiei 

i=l 


Upon  applying  the  operator  eMt,  t>  0,  to  x  we  obtain: 


cMtx  =  £  cxie^itei 
i=  1 


Defining  the  resulting  vector  in  (2.2)  as  x(t),  we  observe  that  the  vector: 


x(t)  = 


|x(t)  | 


can  be  approximated,  for  large  enough  t,  by: 


x(t)  = 


|x(t)  | 


X(t)  =  t  PjCj 


(2.4c 


« 

I 


’  |x(t)| 

and  n  is  a  suitably  chosen  (truncation)  integer.  The  error  involved  in  the  approxima¬ 
tion  (2.4)  is  of  the  order  Je(>^*1_>‘l^J.  Thus,  for  a  given  desired  accuracy  there  are 
values  of  n  and  of  t  which  fulfill  these  requirements.  The  way  t  and  n  are  to  be 
chosen  in  a  practical  algorithm  is  explained  below.  In  the  rest  of  this  section  it  is  as¬ 
sumed  that  (2.4)  is  an  actual  equality  in  order  to  simplify  the  presentation. 

Consider  the  following  vectors: 

vm  =  Mm-1x(t)  (2.5) 

where  1  <  m  <  n+ 1.  Using  (2.4b),  we  obtain: 

vm  =  1*  m  £  n+1  (2.6) 

i=  1 


Consider  now  the  matrix  Wy  defined  as: 

Wy  =  X/*1  1  £  i,  jf£  n  (2.7) 

and  denote  P1ei  =  e,.  The  vectors  ey  1  <  i  <  n  are  n  linearly  independent  vectors 
since  the  e/s  obviously  are.  Eqn.  (2.6)  can  be  rewritten  as: 

vm  =  XWBtkJk  (2.8) 

k=  1 


The  matrix  W  is  a  Vandermonde  matrix  whose  determinant  is  FI  ( A,  -  a.)  and  :s,  by 

»*  j 

assumption,  nonvanishing.  Consequently  the  n  vectors  v,,  are  independent 

and  spanned  by  e,t  e:,  •  •  •  en.  Hence  vn+ j,  which  is  also  spanned  by  Cj,  ■  .c.,,  can 
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be  written  as: 


vn+l  =  Ylvl  +  Y2V2  +  ’  •  •  Ynvn 


Consider  now  the  following  linear  transformation  defined  by: 


Tvm  =  v 


in  ~  vm+ 1 


1  <  m  <  n 


(2.10) 


This  is  a  linear  transformation  from  the  subspace  spanned  by  vlf  •  •  •  vn  or  C\,  ■  •  •  en 
to  itself.  Using  (2.6)  and  the  linearity  of  T: 


Hence: 


Tei=  ZWfii 

1=  1  i=  1 


Tej  = 


(2.11) 


(2.12) 


i.e.  the  X,’s  are  also  eigenvalues  of  the  operator  T  that  acts  in  a  finite  (and  small)  di¬ 


mensional  space.  Let: 


(2.13; 


et=  ZSijvj  1  *  ».j  *  n 

j=i 


be  an  expansion  of  the  ej’s  in  terms  of  the  Vj’s.  Then: 


Tei  =  IS.jTVj 
j=i 


Define  now  the  matrix  D  by: 


Tv,  =  D,.v.  1  <  i,j  <  n 


where  the  (Einstein)  summation  convention  over  repeated  indices  is  assumed.  It  to] 


lows  from  (2.10)  that  D(J  =  5J>1+1  for  1  £  i  <  n-  1  and  from  (2.9;  that  Dn  .  -  y 
D  is  just  a  representation  of  the  T  transformation. 

Upon  substituting  (2.12),  (2.13)  in  the  left  hand  side  of  (2.14),  and  (2.15) 
in  the  right  hand  side  of  (2.14)  we  obtain: 


^■i§ikvk  =  Sij^jkvk 


(2.16) 


Hence: 


^iSik  =  gyDjk 


(2.17) 


i.e.  the  row  gu  (i  fixed)  is  a  left  eigenvector  of  the  matrix  D  corresponding  to  the 
eigenvalue  Using  the  definition  of  the  matrix  D  we  can  now  proceed  to  find  and 
Si.k-  Inspection  of  the  matrix  D  for  small  n  leads  to  the  following  result,  which  can  be 
easily  checked  by  substitution  into  (2.17): 


a.  -  V 

8a  "  U+l-n 
m=  1 


( 2. 1 S ) 


,  y  n+  1- m 
1=  I  A.: 


(2.19' 


Eqn.  (2.19)  can  also  be  rewritten  as: 


I  T.1,-1  -  X»=  0 

m-  l 


(2.20) 


Consequently,  given  the  values  of  yra,  one  can  use  (2.20;  to  compute  the  sr 


;  spe.tr  :m 


and  (2. IS)  and  (2.13)  to  compute  the  eigenvectors.  It  or’'-  remains  :<>  compute  ’he 
values  of  ym.  To  this  end  we  perform  an  orthogonaiication  «>t  the  vectors  v  .  ■  v 


leading  to  the  orthonormal  set  wt  •  •  •  wn  such  that: 


wi=  Xcuvj 
)=  i 


where  the  definition  of  the  matrix  c  is  obvious.  Hence: 


vn+i=  X  ('W"|Ct)wk; 


k=  1 


and  using  (2.21) 


n  n 


vn+ 1  =  I  X('Wwk)Ckmvm 


m=  lk=  1 


(2.21 


'2  22) 


(2.23) 


Hence: 


□ 

y m  ~  X  (vn+l‘wk)ckm 


k=  1 


(2.24) 


Eqns.  (2.18),  (2.20)  and  (2.24)  constitute  the  sought  solution  for  the  eigenvalues  and 
the  eigenvectors. 

All  of  the  above  formalism  can  be  easily  modified  when  the  vectors  vm  are 

defined  by: 

vm  =  e(m‘1)MTv(t)  (2.25 

where  t  is  a  chosen  time  interval.  In  this  case,  in  eqns.  (.2.6)  —  (2.7)  ;s  to  be  re¬ 
placed  by  e''"'1.  This  choice  is  especially  convenient  when  the  matrix  M  actually 
represents  a  differential  operator  (such  as  the  Orr-SommerfelJ  operator).  In  such  a 
case  the  vectors  vm  can  be  obtained  by  solving  x  =  Mx  tor  0  <  i  C  t  w;:  ' 
the  initial  condition,  i.e.  solving  an  initial  value  problem  for  the  differentia’.  "pe:v 


In  this  way  one  avoids  storing  a  large  matrix  M  representing  the  differential  operator. 

Finally,  we  note  that  in  the  above  procedure  the  coefficients  of  e,  for  in¬ 
creasing  values  of  i  are  increasingly  damped  by  exponential  factors.  Moreover,  the 
above  method  should  work,  strictly  speaking,  only  for  cases  of  nondegenerate  spectra. 
A  remedy  to  both  of  these  problems  is  provided  by  a  direct  construction  of  the 
orthogonal  vectors  w;  [cf.  (2.3)  and  (2.21)].  This  is  done  by  defining: 

w  j  =  x(t)  (2.26) 

and: 

i 

Twi=  L  aikwi  +  ai.i-Mwi+i;  is  n  (2.27) 

k=  1 

where:  Twj  =  Mwj  or  Twj  =  eMtWi  corresponding  to  Tej  =  or  Tej  =  e^e^  respec¬ 
tively  (which  are  the  two  alternatives  presented  above).  The  advantage  of  computing 
the  vectors  w  by  the  use  of  (2.27)  is  that  it  enhances  the  weight  of  nonleading  eigen¬ 
vectors  (see  Section  III  for  details). 

The  net  result  of  the  formalism  presented  in  this  section  is  the  construc¬ 
tion  of  n  linearly  independent  vectors  which  are  spanned  (to  a  good  approximation) 
by  the  first  n  eigenvectors  of  M.  Then  either  equation  (2.20)  or  a  representation  of 
the  operator  T  in  the  basis  defined  by  {w,;  1  <  i  <  m]  can  be  used  for  the  computa¬ 
tion  of  the  eigenvalues  and  eigenvectors.  The  resulting  method  is  economical  and 
simple  to  implement. 
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III.  ANALYSIS  OF  THE  METHOD 


» 


* 


The  present  section  is  devoted  to  an  analysis  of  the  formalism  developed  in 
section  II.  More  specifically,  we  investigate  the  error  involved  in  using  a  finite 
number  of  vectors  (or  Wj).  In  what  follows  we  shall  use  the  "exponential"  version 
for  the  operator  T;  Te,  =  e^e,. 


Firstly,  we  wish  to  investigate  the  extent  to  which  the  (right)  eigenvectors, 
(ej,  of  M  are  spanned  by  the  vectors  (wj  or  by  the  vectors  {v;}.  Then  we  propose  to 
estimate  the  error  involved  in  the  computation  of  the  eigenvalues. 


We  wish  to  show  that  the  eigenvectors  e;  can  be  written  as  follows: 


where  R  is  the  rank  of  the  matrix  M,  s  Xj  -  /ij,  t  is  the  time  of  filtering  or  of  ini¬ 
tial  integration  and  the  coefficients  a^  are  of  order  unity.  The  proof  proceeds  by  in¬ 
duction.  By  construction: 


x(t)  = 


R  it 


(3.2) 


where  a,  are  coefficients  [see  eqn.  (2.1)].  Consequently,  wlt  being  a  normalized  vec¬ 
tor  in  the  direction  of  x(t),  can  be  written  as: 


w,  =  N 


R 

lie 

i=  1 


a,e; 


(3.3) 


where  the  right  side  of  eqn.  (3.2)  was  divided  by  e^'1  and  normalized  Nj  is  a  normal¬ 
ization  factor.  Hence: 
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pnnnnnrw  ^  w.  w n  ’J.J J'J  M  «J U ■  Jf  .!>■! >  P I  ij> » l'  ■  tf  ■  1-m  ■  U ■  [ . p .  j i fr  .*■ ■ , „ ^ ^  A ^ „,.v.g.v. 


Ci=  ttV'v>- 

lN  I  al  a!  i=2 


i  3.4) 


a, 


which  shows  that  e,  indeed  has  the  expected  form  (with  an  =  — -  and  a,  = - - 

*  ‘  \'  r4  ‘  *  rj 


N  iai 


for  i  >  2).  Next  we  prove  that  eqn.  (3.2)  is  valid  for  i  =  2.  By  definition: 


w2  =  N2(v2  -  (vj-w/jw,) 


a, 


(3.5) 


where  N2  is  a  normalization  factor.  Hence: 


-^-w2+  (vvw^wj  =  a1e1,<t+t)e1  +  a2eX:(t+T)e2  +  ^  °tlex,(£^x)el  (3.6) 


i=3 


Substituting  ej  from  eqn.  (3.4)  in  the  right  side  of  eqn.  (3.6): 


-~-w2+  (v^wfjWj  a  -j^-eX,(t+T)w1  -  eX,(t+T)eX:,W2 


-  eXl(t+T)^eX,,ta1ei  +  a2eXl(t+T)e2  +  £  a.e^^e,  (3.7) 

i=3  i=3 


By  solving  for  e2  in  eqn.  (3.7)  it  is  easy  to  see  that  formula  (3.1)  is  correct  for  i  =  2. 
Assume  next  that  eqn.  (3.1)  has  been  proven  for  1  £  i  <  m.  Define: 


Jaij  1  S  j  <  i  <  m 
A‘J  =  lo  m  >  j  >  i  >  1 


aue  "  m  >  j  >  i  >  1 
B,i  =  IQ  1  <  j  5  i  <  m 


Both  A  and  B  are  square  matrices  of  rank  m.  Define  also: 


(5.8) 


(3.0) 


a„e 


jc  j  >  m  +  1 ;  1  <  i  <  m 


Q  - 

;i  ~  10  otherwise 


-i: 


By  the  induction  assumption  [and  using  definitions  (3-8,  9,  10)  to  rewrite  eq.  (3.1)! 
we  have: 


e.  =  IA.jwj  +  IVj  +  Icijej 

j* 1  j*  1  J=  1 


for  1  £  i  ^  m.  Hence: 


e,  =  (I  -  B)^1(Akjwj  +  Ckjej) 


where  the  summation  convention  is  assumed.  For  every  integer  r: 


CBr)i>k  =  aUl  aiwil  •  •  •  aif_like 


•  •  •  V,r_,t 


(3.11) 


(3.12) 


(3.13) 


or 


«  e 


(3. 


14' 


Consequently: 

(I  -  B)£l  =  8*  +  (3.15 

where  h^  is  0(1)  and  h^  =  0  for  i  >  k  or  k  >  m  [see  eqn.  (3.9)].  Substituting  eqn. 
(3.15)  in  (3.12)  we  obtain: 

e,  =  AjjWj  +  Cj^j  +  h^e^'A^Wj  +  hlke^“lCkjeJ;  i  <  i  <  m  (3.17 

The  first  and  third  terms  in  the  right  side  of  eqn.  (3.17)  are  linear  combi¬ 


nations  of  {wk;  1  <  k  <  m  }.  The  second  term  is,  using  eqn.  (3.10): 


T"'  A.  ,,t 

X  aue  ■  ej 


The  fourth  term  is: 


R  1  .  R  i  .  R  R  i  | 

X  £  akJe  ikej=  X  (  X  hikak;)e  ”  ej 


(3.19 


k=  i+  1  j=m+l 


j=  m+  1  k=  i+  1 


Consequently,  by  separating  the  terms  containing  {w^l  £  j  <,  i}  and  those  which  are 
superpositions  of  {w^i  +  1  <  j  £  m}  in  eqn.  (3.17)  we  obtain: 


ei=  XYijwj+  X  X  eX*,lhikAitjWj  +  £  rljeXl'teJ 

j=  1  k=i+l  j=i+l  j=  m+ 1 


where  Yy  and  ry  are  0(1)  quantities.  Note  that  we  used  the  fact  that  h^  =  0  for  i  >  k 
or  k  >  m.  It  follows  from  eq.  (3.20)  that 

ei  =  X  Y ijwj  +  X  ®se^wj+  X  rye^ej  (3-^1 

j=  l  j=  i+ 1  j*  m+ 1 

for  1  <  i  <  m,  where  the  definition  of  Gy  is  obvious.  Next  we  use  eqn.  (3.21)  to 
complete  the  induction  process.  Eqn.  (3.21)  itself  follows  from  the  induction  assump¬ 
tion  for  1  <  j  <  m. 

By  its  definition: 


wm+l  =  Nm+1(vm+1-  X(vm+rwi+X) 


where  Nm+1  is  a  normalization  factor.  Consequently: 


vm+  1  =  X 


where  are  0(  1)  numbers.  Using  (2.4b)  and  (2.25): 


m+  1 


Z  qS 


i=l 


.,iwi=  ZP,^”’’*, 
i=  I 


R 

Z  Pi* 

i=  m+  1 


X,(mx)d 


(3.24) 


By  assumption,  the  vectors  1  £  i  £  m  can  be  expressed  by  eq.(3.21).  Substituting 
eqn.  (3.21)  in  eqn.  (3.24)  we  obtain,  after  a  rearrangement  of  terms: 


m+ 1  m  , .  .  R 

ZqSi.i»i=  ZP,<=  ■<”,)  Z 

i=  1  i=  1  j=m+l 


r.je 


R 

V  S 

i=  m+  1 


(3.25) 


where  ^  are  the  new  coefficients  of  the  Wj  terms.  Solving  for  em+1  in  eq.(3.25) 
we  obtain  (using  eq.  (2.4c)  for  3i) : 


m  R 

em+l  =  Z  sm+l,ke 

i=  1  k=  m+2 


4. 


(3.26) 


where  qm+lii  and  sm+1(k  are  0(1)  quantities  (assuming  eX,jt  is  0(1)).  This  completes 
the  induction.  Consequently  eqn.  (3.1)  and  eqn.  (3.21)  are  correct  for  all  i  >  1. 


It  follows  from  eq.  (3.21)  that  the  error  involved  in  the  assumption  that  e, 
is  a  combination  of  Wj,  w2,  •  •  •  ,  wm  is  0(e^B*l,,t),  which  means  the  choice  of  the  size 
of  the  subspace  of  (w^  1  ^  i  <,  m)  should  be  such  as  to  have  !Re(Xm*  1  -  Xj)t 1  >  >  1 
for  1  <  i  <  r  if  r  correct  eigenvectors  are  wanted.  Notice  that  no  gap  in  the  spectrum 
is  necessary  for  this  estimate  or  the  proposed  method  to  be  valid.  Some  of  the  eigen¬ 
vectors  (i.e.  those  for  which  ;Re(A.m+l-  \j)t|  is  not  large  enough)  will  not  be  well  ap¬ 
proximated  by  the  procedure.  In  such  a  case  an  appropriate  increase  ot  the  value  ot  m 
will  ensure  a  good  approximation  for  er  A  similar  statement  will  be  shovsn  below  to 
be  true  for  the  corresponding  eigenvalues.  Before  we  do  that,  we  present  a  second 
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result  approximating  the  error  involved  in  expressing  the  eigenvectors  e,  in  terms  of 
the  vm’s  [see  eqns.  (2.8),  (2.9)]. 


We  wish  to  show  that  for  each  1  <  i  <  R  there  is  a 
{vp,  1  £  i  S  R},  which  we  denote  by  Uj  and  which  satisfies: 


ui  =  ei+  £  Kue>Lj’tej 
j=1+i 


where  Ky  are  0(1).  As  before  this  statement  may  be  proven 
duction.  Here  we  shall  only  sketch  a  proof.  For  i  =  1  define 


ui  = 


“i=  i=i 


Hence: 


R  a;  ,  . 

“i  =  *.+  s 

j=2  U1 


Next  define  as  vk  divided  by  a1e*'l(t+^-  l^: 


vk<»»  «,+  2— LeX’,<,*“"1)” 

j=  2  al 


Hence: 


R  a 


'k  -«i  =  I—  e  “ 

j=  2  al 


J  ^>-,1  1)T 


-l)ej 


a-.  , 


Dividing  \v(!)-u.  by  — we  obtain: 

K  al 


linear  combination  of 

(3.27) 

by  straightforward  in- 

(3.28) 

(3.29) 

(3.30) 

(3.31) 


6- 


vk2)  =  V—  e*-,,‘cei',<k_  l,T—  1  )Cj 


Next  define: 


u2  -  v|2V(e^2lt  -  1) 


(3 


) 


or: 


u2 


e2  + 


R  otj  e^(e^  -  1) 
j=3  «2  e^*x  -  1  Cj 


A  similar  procedure  of  Gauss  elimination  steps  leads  to: 


(3.33) 


u3  = 


R  a 

e3+  £-<-e 
j=4  a3 


e2Xi,T-  1  eXil-c-  1 

e2X2,-c_  1  e^i_  {  ^ 

e2*,lX-  1  e^lX- 1 

e2)l2,t-  1  eXjlT-  1 


(3.34) 


The  continuation  of  this  process  obviously  leads  to  the  desired  equation 
(3.27).  It  is  now  easy  to  see  the.ej  can  be  expressed  in  terms  of  (v^  •  •  •  vm)  with  an 
error  of  0(eX'ffl*’,it).  Notice  that  when  x  itself  is  large  the  coefficient  of  ej  in  eqn.  (3.34) 
tends  to  zero,  which  shows  that  one  may  obtain  adequate  approximations  for  large 
values  of  i  and  (even)  short  values  of  t. 

Finally  we  turn  to  estimating  the  error  in  the  eigenvalues,  when  computed 
using  a  finite  number  (n)  of  vectors.  To  this  end  define  a  matrix  U  by: 

U.j  =  WjTWj  1  <  i,j  <  n  < 3  35 

U  is  the  reduction  (or  projection)  of  the  operator  T  to  the  finite  d::r. en- 
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sional  subspace  spanned  by  the  v’s  (or  w’s).  Let  <{>  be  an  eigenvector  of  L',  with  a 
corresponding  eigenvalue  A: 


=  A4>i 
j=i 

Let: 


ei  =  L  bijwj  1  ^  i  ^  R 
j=i 


Then,  since  Tej  =  e  1  eif  we  have: 

R  XrR 

t  t>ijTwj  =  e*'  £  byWj 

j=i  j=i 


Hence: 


Twi=  Y, 
j.*=i 


kiwj 


and,  using  eqn.  (3.35): 

Uu  = 

k=  1 


Hence,  from  eqns.  (3.36)  and  (3.40): 


A*; 

j=  lk=  1 


( 3.36 


(3.37 


(3.38 


'3.3$ 


( 3.4C 


(3.41 


Define: 


rmrmi  ^ F  JVT  *A.W  v  ^  KT  *.^  *Lr"  HJ*  H."  HJ*  lor  *."  HWV  V' 


VWWV  V 


I 


I 


Hence 


X  X  bribik'e^'Vic  =  Ayr 
i=  lk=  1 


(3.43) 


Since  Xbrib  &  =  8rk,  we  have: 
i=  1 


R  R 

=X,V,  +  £  I  bribi'e**'^  A% 

k=  1  i=m+  1 


(3.44) 


By  comparing  eqn.  (3.37)  and  eqn.  (3.1)  or  (3.21)  we  observe  that: 


b;i  = 


lij  j  £  i 

vv  j  >  ‘ 


(3.45) 


where  the  definition  of  a^  is  obvious.  Using  eqns.  (3.8),  (3.9)  with  m  =  R  and  a,;;  re¬ 


placed  by  a;;  in  eq.  (3.9)  we  see  that: 


b  =  A  +  B 


Hence: 


(3.46) 


b"1  =  A-1  -  A^BA"1  +  A-1BA-1BA-1 


(3.47 ) 


Since  Akl'  =  0  for  k  <  i,  by  definition,  we  have  to  estimate  only  terms  containing  ele¬ 
ments  of  the  matrix  B  in  (3.47).  In  this  case: 


(A_lBA_I)lk-  A~ ’Bj^A"^ 


( 3.48 ) 


or: 


-19- 


!  ? .  4 ' )  1 


(A-1BA-')lk=  I  A"1  V  a„1:e  *  A"^ 

i,<  i  ii<  i: 

The  largest  contribution  of  e  Vl  which  is  possible  under  the  constraints 
k  <  i2  >  ij  ^  i  is  for  i2  =  k  and  ij  =  i.  A  similar  analysis  holds  for  all  terms  in  the 
sum  eq.  (3.47).  Hence: 

1  =  eXk,tdiJC  for  i  <  k  (3.50) 

where  are  0(1)  numbers.  When  i  >  k,  A^1  is  0(1)  itself  and  we  denote  it  by  d*. 
Hence,  from  (3.45): 

R  i  R  R 

eX,Vr+  £  £  Mik^'Vk  +  £  £bndA\=AVr  (3.51) 

i=m+lk=l  i=m+ lk=i+ 1 

Using  (3.45)  again  (to  substitute  for  bri): 

+  £  £  Tne^d^e^Vk  + 

i=  m+  1  k=  1 

R  R 

£  £  arieX*tdikeX‘Vk'Vk  =  Ayr  (3.52) 

i=m+l  k=i+l 


For  |Re(Xm+1  -  Xr) 1 1  >  >  1  the  leading  order  terms  in  the  double  sums  of  eqn. 


(3.52)  are  eXmr'J\m+l  £  dm+  ltke*k>k  and  c'l,*8JX.m+1da+lim+2e'B,:1\|/a*I,  respec- 

k=  1 


m*  1 


-  . . .  _ 1  j  -  ^  - 


tively.  Hence: 


(e^1  -  A)vr  +  e^*1-'1 


m+  1 

,m-<- 1(  £  ^ 
k=  1 


m+l.ke  Hfk 


+  e' 


n*  i 


1  ,m+ . 


Mf- 


.)  =  <J<  3.5! 


To  lowest  order  (i.e.  neglecting  e^m*l,‘)  an  eigensolu tion  satislies:  A  =  e"“  and 
i/j  =  5jr.  The  next  perturbative  correction  yields: 
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which  indeed  shows  that  for  a  large  enough  filtration  time,  the  error  in  the  eigen¬ 
values  is  exponentially  small  provided  ReAm+1-t  is  well  separated  from  (or  much 
smaller  than)  ReXr-t.  Thus  no  gap  in  the  spectrum  is  necessary  to  obtain  excellent 
values  for  the  leading  eigenvalues.  Obviously  an  increase  in  m  will  make  the  term 
e^-u1  as  smai]  as  needed. 


it 


IV.  THE  ALGORITHMS 


In  this  section  we  present  detailed  algorithms  for  the  computation  o 
/ij,  •  ■  •  ,X0  an^  ,en.  The  algorithm  is  described  in  steps: 

(1)  Choose  an  initial  vector  x0.  It  does  not  matter  if  x0  is  independent  of  Oj 
or  other  relevant  eigenvectors.  They  will  be  introduced  in  step  (2)  by 
round-off  errors  during  the  filtration  process. 

(2)  Solve  the  equation  x  =  Mx,  with  x0  as  initial  condition,  up  to  a  time  t. 

x  ( t) 

Normalize  the  resulting  x(t):  x,  =  — - .  Use  now  x,  as  an  initial 

1  Ix(t)|  1 

condition  and  solve  for  a  time  8,  obtaining  x^Q).  Repeat  this  process  r 
times  to  obtain: 


x(r9) 
|x(r0)  | 


(4.1 


The  reason  one  normalizes  after  each  time  t  is  to  avoid  deal¬ 
ing  with  large  numbers  (since  |x(t)  |  «  e^’1  for  large  times).  The  choice 
of  r  is  explained  below.  Since  it  is  not  clear  apriori  whether  x0  is  in¬ 
dependent  of  ej,e2  etc.,  it  is  advisable  to  use  low  accuracy  computation 
in  the  first  few  iterations.  In  this  way  the  weight  of  e^ei,  etc.,  will  be 
amplified. 

(3)  Compute  vk  =  Mk_1xn,  1  <  k  ^  m+1.  Orthonormalize 
Vj,  v:,  •  ■  •  vm  [using  Householder  transformations  (Wilkinson,  1965), 
for  example!.  The  resulting  orthonormal  vectors  are  denoted  by 
Wj ,  \v;,  •  •  •  Wm.  Following  each  step  of  the  onhonormalizution  pro- 


'S 


cedure  use  the  test  of  step  (4)  to  make  sure  v,.  •  ■  •  v,r  are  independent. 

i 

The  matrix  c,  defined  by  \v;  =  Yc^v,  that  results  from  the  orthonormal- 

j  =  1 

ization  procedure  should  be  kept. 

(4)  Test  whether  vm+1  is  spanned  by  v]t  •  •  •  vm.  To  do  this  compute: 

HE  !!  =  |vB+1-  |vm+1|  (  4.2} 

i=  l 

If  the  error  E  is  larger  than  desired,  increase  m.  If  this  process  results  in 
too  large  a  value  for  m,  go  back  to  step  (2)  and  do  several  additional 
iterations.  Then  repeat  steps  (3)  and  (4)  until  E  is  smaller  than  the 
desired  accuracy. 

(5)  Once  the  matrix  c  and  the  vectors  vk  and  wk  are  known,  use  (2.24)  to 
find  the  yk’s.  Alternatively  solve  the  least  square  problem  minimizing 

the  expression  ||vm+1  -  2,Yivjll  •  Standard  least  square  routines  may  be 

i=l 

used  to  solve  for  the  y;’s  directly. 

(6)  Use  (2.20)  to  find  the  spectrum. 

(7)  Use  (2.18)  to  find  the  eigenvectors. 

It  should  be  stressed  that  in  step  (4)  the  test  should  be  performed  for 
m  =  2,  3,  etc.,  up  to  a  desired  m  so  as  to  make  sure  that  vt,  ■  •  •  va  are  indeed  in¬ 
dependent.  It  may  happen  that  the  test  in  step  (4)  results  in  a  value  of  m  which  is 
smaller  than  its  desired  value.  In  this  case,  either  use  a  shorter  integration  time  t  or 
the  procedure  described  below.  If  spurious  eigenvalues  appear,  they  will  be  m- 
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dependent.  A  comparison  of  the  eigenvalues  for  different  \alues  of  m  enables  the 
identification  of  the  nonspurious  eigenvalues. 


The  reason  the  procedure  described  above  may  converge  for  a  vaiue  of  m 
(say  mj)  which  is  smaller  than  the  desired  (number  of  eigenvalues)  m  can  be  a  gap  in 
the  spectrum  at  m]t  namely: 

Re^!  >  ReX2  •  •  ■  ReXm|  >>  Re\mi+1  >  •  •  •  (4.3) 

In  this  case,  Xj,  ‘  '  '  an<^  corresponding  eigenvectors  are  obtained  to  an  ex¬ 
ponential  accuracy  which  is  demonstrated  in  the  next  section.  Therefore  it  is  possible 
to  obtain  the  mi  left  eigenvectors  e,L  corresponding  to  ei(  1  £  i  £  m ^  These  vectors 
satisfy: 

e,L  e,  =  5if  (4,4) 

Expanding  ej  in  the  orthonormal  set  w  (see  eqn.  (2.21)): 

(e,L-wk)(wk+ej)  =  Sjj  (4.5) 

Thus  (e1L-wk)  is  the  inverse  of  the  (known)  matrix  A  (of  rank  m^: 

Qkj  =  (wk+ej)  (4-6> 

Consequently: 

e,L=  iV'")'  (4'’ 

!■  1 

Consider  now  a  vector  x0  that  is  chosen  to  be  orthogonal  to  e,L,  1  £  i  ^  m.  Now 
perform  step  (2)  with  the  modification  that  after  each  time  0,  xr(t)  is  orthogonalized 
to  Wj,  •  •  •  wm  The  resulting  vector,  after  a  time  t.  is: 


The  difference  [ej- (ejWj^Wj]  is  exponentially  small  for  the  exact  e/s.  So  is  ( e, w;+ ) 

for  i  >  Thus,  for  t  not  too  large:  y(t)  «  £  e^cq  and  the  regular  procedure  can 

i>  m) 

be  applied  to  y(t)  to  find  Xm)>  etc.  If  Xm)  is  close  to  Xmi+1  a  relatively  large 

time  t  in  (4.10)  may  be  necessary.  In  this  case  a  better  accuracy  for  ej,  1  £  i  <  m, 
should  be  obtained  first.  Methods  to  do  so  will  be  given  in  a  future  publication. 

Another  version  of  the  suggested  algorithm  which  is  highly  suitable  when 
no  gap  in  the  spectrum  exists  is  based  on  eqns.  (2.25)  —  (2.26).  Step  (2)  involves 
then  an  orthogonalization  after  each  integration  for  a  time  of  x.  The  projection  of  the 
operator  T  (in  the  subspace  spanned  by  the  w’s)  can  then  be  expressed  as  shown  in 
Section  III  [eqn.  (3.35)].  The  small  resulting  mxm  matrix  can  be  diagonalized  by 
LR/QR  methods  (Wilkinson,  1965)  (we  have  used  EIGCG1.  an  EISPACK  routine! 
yielding  exponentially  accurate  results. 

We  summarize  this  section  by  mentioning  that  the  above  method  can  be 
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V.  IMPLEMENTATION  OF  THE  METHOD:  THE  ORR-SOMMER FEED  EOLA¬ 
TION. 

In  this  section  we  present  an  application  of  the  methods  described  above  to 
a  problem  in  hydrodynamical  stability":  the  spectrum  of  the  Orr-Sommerfe!d  equation 
for  channel  flow.  This  equation  has  an  infinite  number  of  degrees  of  freedom  and 
there  is  no  "gap"  in  its  spectrum.  Thus  it  is  a  good  test  case  for  our  methods.  More¬ 
over,  the  existence11  of  previously  computed,  highly  accurate  values  of  the  spectrum 
enables  us  to  perform  a  comparison  of  our  results  with  the  known  spectrum. 

The  Orr-Sommerfeld  equation  reads: 


(U-c)(D2-cr)y-U'Y  =  — (D2-cr)2v  ;5.1) 

iaR 

V(x.y.t)  =  V(y)exp(ia(x-ct)J 

where  UCy)  for  -  1<  1,  is  the  basic  velocity  profile,  v(x,v,t)  is  the  perturbation 

d 

streamfunction,  R  the  Revnolds  number  and  D  the  cross-stream  derivative  — — .  The 

dy 

real  number  a  represents  the  wavenumber  of  the  streamwise  periodic  perturbation 
and  c  is  the  sought  complex  eigenvalue.  The  real  part  of  the  eigenvalue  c,  is  the  phase 
speed  of  the  perturbation.  The  growth  rate  of  the  perturbation  is  exp:  etc. t)  where 
c,  =  Im  c.  Thus,  a  perturbation  v  is  stable  if  c,<  0  and  unstable  if  c,>  0.  The  boun¬ 
dary  conditions  on  eq.(5-l)  are: 

V<±  1)  =  Dy(±  1)  =  0 

The  above  is  an  eigenvalue  problem  It  is  transforme  a  time 
dent  problem  {  from  winch  it  actually  derived)  by  redefining  the  •  .  e  *> . 


function  y  as: 


V(x,y,t)  =  y(y,t)e10U 

even  though  the  equation  is  separable  in  time.  We  thus  solve  the  following  initial 
value  problem: 

3t(D2-or)\j/+  iaU(D2-ot2)y-  iaU"y  =  R_1(D2-a2)2y 


y(±  1)  =  Dy(±  1)  =  0 

( 5.2a.b) 

It  is  convenient  to  define: 

C  =  (D2-or)y 

(5.3) 

^  =  p+a 

(5.4a) 

V  =  <t>  +  iX 

(5.4b) 

to  obtain  the  coupled  real  system: 

3tp-aUX+aU"x  =  R_1(D2-a2)p 

(5.5a) 

9tX+aUp-aU"0  =  R-^D^crA 

(5.5b) 

(D2-a2)0  =  p 

(5.5c  i 

(D  2-cr)X  =  ^ 

( 5.5d ) 

with  boundary  conditions: 

0(±  1)  -  X(-  I)  ~  D0(±  1)  =  Dxf-  1  '>  =  V 

i  5  fi 

The  real  system  may  be  formally  written  as  y  =  My,  M  being  an  integro-ditferentul 
operator. 

These  equations  have  been  solved  by  the  Chebycheff  pseudo-spectral  tech¬ 
nique  (see  the  appendix  for  details  ).  An  initial  filtering  time  Tf  was  used  to  produce  a 
function  y(y,t).  which  is  spanned  by  a  relatively  small  number  of  eigenvectors 
corresponding  to  leading  eigenvalues  Re(-iac)  or  the  fastest  growing  modes,  as  has 
been  explained  previously. 

In  all  computations  results  have  been  presented  for  the  parameter  values 
a  =  1.0  and  R  =  10,000.  The  base  velocity  profile  U  is  (1-  y2). 


The  results  of  our  computations  using  the  methods  developed  are  present¬ 
ed  in  five  tables.  A  sixth  table  of  eigenvalues  from  refill  is  taken  as  standard  and 
used  for  comparing  the  accuracy  of  our  results  In  all  cases  we  obtain  several  leading 
eigenvalues  to  a  very  high  accuracy. 

In  our  computations  we  vary  the  filtration  time  Tf,  the  sampling  interval  t, 
the  accuracy  of  the  time  integration  as  expressed  by  the  discrete  time-step  size  during 
filtration  5tf  and  during  the  sampling  diz  in  order  to  elucidate  their  effect  on  the  accu¬ 
racy  of  the  computed  spectra  and  the  possible  generation  of  spurious  modes.  We  have 
also  investigated  the  influence  of  the  initial  streamfunction  y(y,0l  and  of  the  number 
of  eigenvalues  determined  on  the  accuracy  of  the  results. 


Three  methods  have  been  used  for  obtaining  the  eigenvalues,  as  exp. amen 
below.  Table  I  presents  results  of  2  runs  with  a  sampling  time/delay  time  of  la  and 
an  initial  filtration  of  50.  The  number  of  correct  eigenvalues  obtained  using  eq.i  2-2"; 


is  5  when  the  number  of  vectors  is  S  or  16.  The  reason  for  this  is  the  hut  that  ’''e 
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filtration  time  is  too  large  to  produce  more  eigenvalues  accurately.  For  the  case  n~  lo 
we  also  obtain  two  dominant  spurious  modes.  This  is  due  to  the  sensitivity  of  the 
roots  of  the  characteristic  polynomial  to  small  inaccuracies  in  the  coefficients.  The 
dominant  spurious  eigenvalues  may  be  easily  identified  by  their  appearance  as  n,  the 

number  of  vectors  is  increased.  The  initial  condition  for  table  I  is  (l-y)"-  V  y,, 

i=  1 

where  y)  are  eigenvectors  corresponding  to  the  fastest  growing  8  eigenmodes. 

Table  II  was  obtained  with  an  initial  condition  of  y  =  U2(l+y).  The  factor 
(1+y)  is  introduced  to  ensure  that  the  initial  condition  has  both  an  even  and  an  odd 
part  as  U“  is  even  in  y.  As  the  filtration  time  is  increased  from  Tf  =  50  to  Tf  =  75  the 
number  of  accurate  eigenvalues  rises  to  5.  The  longer  filtration  time  allows  for  a 
better  damping  of  the  decaying  modes  not  of  interest.  It  is  possible  that  the  cruder 
time  step  in  the  second  case,  8tf  =  0.075  as  opposed  to  5tf  =  0.05  may  introduce 
modes  independent  of  the  initial  condition  by  way  of  numerical  noise. 

In  Table  III  the  vectors  vm  were  generated  from  the  initial  condition 
U2(la-0.71y)  after  a  filtration  of  Tf  =  75  and  a  delay  interval  of  z  =  20.  How'ever,  in¬ 
stead  of  using  eq.(2-20)  for  finding  the  eigenvalues,  we  have  employed  the  LR  algo¬ 
rithm  \s  implemented  in  the  EISPACK  routine  EIGCG1  to  calculate  the  eigenvalues 
of  the  matrix  directly.  Four  eigenvalues  are  obtained  accurately  both  tor  n=  8  and 
n=  16  vectors. 

Tables  IV  and  V  were  produced  using  the  method  of  orthogonalizauon 
described  at  the  end  of  section  II  (  see  eq. (2-25, 26)  ).  The  eigenvalues  of  the  reduced 
matrix  obtained  were  again  computed  using  EIGCGl.  In  table  IV  an  increase  of  1, 
from  50  to  150  leads  to  an  increase  in  the  number  of  accurate  eigenvalues  from  3  to  n 
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for  n  =  8  vectors.  The  improved  accuracy  using  this  last  method  is  obtained,  as  ex¬ 
plained  in  section  III,  by  the  amplification  of  the  weight  of  the  non-leading  eigenvec¬ 
tors  in  the  orthogonalization  procedure.  Indeed,  as  seen  in  table  V,  the  number  of  ac¬ 
curate  eigenvalues  when  n=  32  vectors  is  used  is  13.  The  method  of  orthogonalization 
produces  12  accurate  eigenvalues  for  n=32  even  when  no  filtration  is  invoked.  In  all 
cases  we  have  obtained  excellent  accuracy  for  the  leading  eigenvalue.  The  subdom¬ 
inant  eigenvalues  corresponding  to  anti-symmetric  modes  tend  to  be  somewhat  more 
accurate  than  eigenvalues  corresponding  to  symmetric  modes. 


VI.  SUMMARY  ANI)  CONCLUSIONS 


We  have  shown  how  one  can  obtain  leading  eigenvalues  and  eigenvectors 
for  large  asymmetric  matrices  using  a  relatively  simple  and  economical  numerical 
scheme.  We  have  also  shown  how  such  a  method  can  be  applied  when  the  matrix  to 
be  diagonalized  is  a  differential  operator.  In  the  latter  case  our  method  does  not  re¬ 
quire  storing  of  an  effective  matrix,  which  represents  the  differential  operator  on  a 
complete  basis  of  expansion  functions. 

Three  variants  of  the  method  were  tested:  a)  A  direct  application  of  the 
formalism  presented  in  section  II.  b)  A  construction  of  a  set  of  vectors  as  explained  in 
section  II  and  the  associated  reduced  matrix,  followed  by  a  diagonalization  using  stan¬ 
dard  eigenvalue  routines  for  general  matrices,  c)  Direct  construction  of  a  set  of 
orthogonal  vectors  as  at  the  end  of  section  II  followed  by  a  diagonalization  of  the  re¬ 
duced  matrix. 

The  third  variant  seems  to  be  the  superior  method.  The  method  was 
found  to  be  very  robust  and  did  not  require  fine  tuning  to  improve  the  accuracy  of  the 
calculated  eigenvalues.  Nor  does  it  produce  spurious  dominant  eigenvalues.  The  accu¬ 
racy  of  the  leading  eigenvalues  may  be  further  increased  by  considering  larger  reduced 
matrices,  that  is  by  increasing  the  number  of  orthogonal  eigenvectors.  The  errors  in¬ 
volved  in  the  proposed  algorithm  are  analyzed  in  detail  in  section  III  and  are  in  good 
agreement  with  our  numerical  results. 

When  a  large  matrix,  of  rank  R.  is  considered  the  number  of  operations 
necessary  to  find  its  eigenvalues,  using  standard  methods,  scales  like  R3.  In  our 
method,  the  number  of  operations  is  proportional  to  the  rank  R,  the  number  ot  vec- 


s 


tors  used,  m,  and  the  time  of  filtration  tf,  i.e.  it  is  Rntf.  The  time  of  filtration  depends 
on  the  spectrum,  as  explained  in  section  III.  If  one  wishes  to  obtain  k  correct  eigen¬ 
values  one  needs: 


t: 

I 

o 

V* 


l- 


Y 


R 


J  <<;  j 

Thus  the  value  of  tf  is  of  the  order  of  — - - — - — . 

I^m+  1  —  I 

Future  applications  of  this  method  could  include  complex  and  non- 
Newtonian  hydrodynamical  stability  problems,  lattice  eigenvalue  problems  and  other 
systems  leading  to  large  asymmetric  matrices  for  which  the  dominant  eigenvalues  are 
of  interest  and  standard  methods  are  either  slowly  converging  or  the  memory  require¬ 
ments  are  prohibitively  large.  We  believe  that  this  method  coupled  to  acceleration 
techniques  may  enable  one  to  tackle  many  interesting  eigenvalue  problems. 
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APPENDIX 


In  this  Appendix  we  present  some  details  concerning  the  numerical  solu¬ 
tion  of  eqn.  (5.4).  We  have  employed  a  pseudospectral  technique  based  on  the  expan¬ 
sion  of  V  in  Chebyshev  polynomials,  which  gives  a  good  resolution  of  boundary  and 
critical  layers.  A  third  order  Adams-Bashforth  time-stepping  scheme  was  used  for  an 
explicit  evaluation  of  the  advection  (variable  coefficient)  terms  and  a  second  order 
Crank-Nicholson  scheme  for  the  linear  diffusion  terms.  The  time  discretized  equa¬ 
tions  read: 


,(n+l)_  n(n) 


- — £_!  _  aU[2ix(°)  -  JAx(n_I)  +  — x(n-2)] 

At  12  12  12 


+  aU"[^|x<0>-  -jfx""'’*  Yrx'”‘2,]=  ~*(D2-  a2)(p(n- 11  -  p(n))  ( A.l) 


♦  au[— p(n)  -  d4p<-»  + 

At  12^  12  12 


-  aU"[21(^(n)  _  il^(n-l)  +  -A.0(n-2)j  _  _L(D2-  a-)[X(0+1)  +  X(n)](A.2) 


Eqns.  (A.l)  and  (A. 2)  can  be  rearranged  as  follows: 


(D2_  a2_  iIL)p(n+l)  _  pjn.n-l.n-2) 


(A. 3) 


(D2_  a2_  2JL)\ln*»  =  pfn.n-2.n-2) 

At 


Fj  and  F->  being  found  from  (A.l)  and  (A. 2).  The  Junctions  o  and  x  satis!>  isce 


■  .*•  Vv*  s ^  O  »>  ■  V  .  *  . 


V  V  v.  /_  -  . 
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eqns.  (5.3),  (5.5)]: 


(D2  -  a2)p(D+1)  =  p(n+1) 


(D2  -  a2)x(n+1)  =  X(n+1) 


(A. 4) 


Since  the  boundary  conditions  are  on  \\f  =  <5  +  i%  alone,  we  have  used  the  following 
Green  function  technique  to  satisfy  these  conditions.  First,  we  find  solutions  to  (A. 3) 
satisfying  p(n+1^(±  1)  =  X^n+1^(±  1)  =  0.  We  call  these  solutions  p^Sm^  and  Ah("* 
These  solutions  are  then  substituted  in  the  right  hand  side  of  eqns.  (A. 4)  and  solved 
using  the  (given)  boundary  conditions  4> (dt  1)  =  %(±  1)  =  0.  We  call  these  solution 
^hom  anc*  Xhoin  respectively.  They  do  not  necessarily  satisfy  the  Neumann  boundary 
conditions  D$(±  1)  =  D%(±  1)  =  0.  Next  we  solve  the  homogeneous  equation: 


(D2  -  a2  -  —)p  =  0;  (D2  -  a2  -  —  )X  =  0 
At  K  At 


(A. 5) 


using  boundary  conditions  p(  1)  =  1,  p(-  1)^0  and  1)  =  1,  X(-  1)  =  0  respectively. 
These  solutions  are  called  p+  and  X+.  Similarly  we  find  p_  and  which  satisfy  (A. 5) 
with  p_(l)  =  0,  p_(0)  =  1  and  X_(l)  =  0,  X_(0)  =  1.  Subsequently  eqns.  (A. 4)  are 
solved  with  p+,  p_  and  X+,  on  the  right  side  yielding  <>±  and  x±  respectively  (with 
boundary  conditions  p±(±  1)  =  X±  (-  1)  =  0  as  required).  The  general  solution  for 
^(n+i)  ancj  ^(n+i)  can  now  written: 


❖  hSrn0  +  a+<|»+  +  a.p. 


=  Xhoml)  +  KX+  +  b_X- 


( A. 6) 


The  constants  a*  and  b*  are  determined  by  imposing  the  Neumann  boundary  corci- 


tions  D0(n+1)(±  1)  =  Dx(n+1)(i  1)  =  0  thus  yielding  a  solution  O'14-1  +  ixn+ 1  which 
satisfies  all  four  boundary  conditions.  The  solutions  $±  and  x±  need  only  be  com¬ 
puted  once  in  a  preprocessing  step  thus  necessitating  only  two  Poisson  solvers  per 
time  step.  The  above  leads  to  a  very  efficient  time  integration  for  the  initial  value 
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Table  I 


Tf  =  50  1  =  15  8tf  =  0.05  5tt  =  0.005  Tf  =  50  x  =  15  5tf  =  0.05  5u  =  0.005 


Identified 

Mode  Number 

Eigenvalue  q 

Identified 

Mode  Number 

Eigenvalue  c, 

1 

3.73967060-03 

1 

3.73967060-03 

2 

-3.51600679-02 

♦ 

4.5S986340-04 

3 

-3.52076045-02 

* 

3.40898714-04 

4 

-5.08987828-02 

* 

-3.02658282-03 

5 

-6.31504249-02 

* 

-1.49838951-02 

* 

-8.72158010-02 

★ 

-2.36492604-02 

* 

-1.24740808-01 

* 

-2.93611665-02 

* 

-1.32275383-01 

2 

-3.51423731-02 

3 

-3.52793842-02 

4 

-5.08971644-02 

★ 

-1.09428509-01 

* 

-1.23518682-01 

* 

-1.48108621-01 

* 

-1.84825899-01 

Tabic  II 


Tf=  5 


Identified 
Mode  Number 


1 

2 
4 


Tr  =  75  z  =  20  5t,  =  0.075  St 


Eigenvalue  cf 


3.73967054-03 

-3.51921344-02 

-5.09056673-02 

-5.52082215-02 

-7.13017236-02 

-1.32670178-01 

-1.47519731-01 

-1.64393611-01 


Identified 

Mode  Number 

Eigenvalue  c, 

3.73967055-03 
-3.51630043-02 
-3.67725906-02 
5.08973087-02 
6.18479421-0 
8.98244690-0 
1.25561448-01 
1.90016840-01 


r  i  r  i 


Table  III 


Tf  =  75  x  =  20  5tf  =  0.075  StT  =  0.005 


Tf  =  50  t  =  20  5tf  =  0.05  5tT  =  0.005 


Identified 

Mode  Number 

Eigenvalue  q 

Identified 
Mode  Number 


3.7396706  IE-03 
-3.513133 10E-02 
-3.77382799E-02 
-5.08939102E-02 
-6.08654660E-02 
-8.94I64908E-02 
-1.23537573E-01 
-I.83295676E-0I 


Eigenvalue  c, 


5.042I7689E-03 
3.73967060E-03 
-9.48657987E-04 
-1. 15472857E-02 
-1.55870014E-02 
-2.303041 83E-02 
-3.40328560E-02 
-3.51905276E-02 
-5.09188577E-02 
-5.121 19074E-02 
-6.27014559E-02 
-9.9343823  IE-02 
-1.1 1875376E-01 
-1.11 94141 6E-0 1 
-1. 17783465E-01 
-1.70313262E-01 


Table  IV 


T,  =  50  x  =  20  Tf  =  75  x  =  20  Tf  =  150  T  =  20 

5tf  =  0.05  5tx  =  0.005  8tf  =  0.075  6u  =  0.005  ot’f  =  0.075  5t.  =  0.005 


Identified 

Mode  Number 

Eigenvalue  c, 

Iden  titled 

Mode  Number 

Eigenvalue  c, 

Identified 

Mode  Number 

Eigenvalue  c, 

1 

3.73967060E-03 

1 

3.73967060E-03 

1 

3.73967060E-03 

2 

-3.519321 16E-02 

'J 

-3.51 137133E-02 

-3.51723718E-Q2 

4 

-5.0905337  IE-02 

3 

-3.74768214E-02 

3 

-3.51870632E-02 ! 

* 

-5.48854885E-02 

4 

-5.08944600E-02 

4 

-5.08987622E-02 

* 

-7.369 1573  IE-02 

5 

-6.07586848E-02 

5 

-6.30009355E-02 

* 

-1.38156654E-01 

7 

-9.17031242E-02 

* 

- 1 .37  863 1 62E-0 1  ! 

♦ 

-1.413711 17E-0I 

* 

-1.23567838E-01 

* 

-1.43432356E-01 

* 

-1.58274564E-01 

* 

-1.82561894E-01 

* 

-1.747288 12E-01 

i 

Table  V 


Tf  =  75  t  =  20  Tf  =  50  t  =  20  Tf  =  0  r  -  20 

5tf  =  0.075  5tT  =  0.005  5tf  =  0.05  5tt  =  0.005  5tf  =  5u  =  0.005 


Identified 

Mode  Number 

Eigenvalue  c. 

Identified 

Mode  Number 

Eigenvalue  c, 

Identified 

Mode  Number 

Eigenvalue  c. 

1 

3.73982608E-03 

1 

3.73967060E-03 

1 

3.73967061 E-03 

2 

-3.51 127796E-02 

2 

-3.51672225E-02 

2 

-3.51672225E-02 

3 

-3.51320836E-02 

3 

-3.51865287E-02 

3 

-3.51865287E-02 

4 

-5.08984842E-02 

4 

-5.0898730SE-02 

4 

-5.08987305E-02 

5 

-6.31517698E-02 

5 

-6.32014442E-02 

5 

-6.32014442E-02 

6 

-6.32018420E-02 

6 

-6.3251521 9E-02 

6 

-6.32515219E-02 

7 

-9.1 1772541E-02 

7 

-9.12220344E-02 

7 

-9.12220344E-02 

8 

-9.12687323E-02 

8 

-9. 1 3 1 34832E-02 

8 

-9.13134832E-02 

9 

-1.19161983E-01 

9 

-1.1922328SE-01 

9 

-1.19223290E-01 

10 

-1.19342042E-01 

10 

-1.19380224E-01 

10 

-1.19380226E-01 

* 

-1.21 170936E-01 

11 

- 1 .245008 10E-01 

11 

-1.24500810E-01 

11 

-1.24S00425E-01 

12 

-1.38224635E-01 

12 

-1.38224635E-01 

12 

-1.38223897E-01 

13 

-1.45450671E-01 

13 

-1.45453318E-01 

* 

-1.44467829E-01 

* 

-1.5106341  IE-01 

* 

-1.5116 1943  E-0 1 

13 

-I.45609239E-01 

* 

-1.58158173E-01 

• 

-  1.58054756E-01 

14 

-1.46448075E-01 

♦ 

-1.60717285E-01 

* 

-1.6071 1699E-01 

15 

-1.75078430E-01 

* 

-1.61 124155E-01 

* 

-1.61259736E-01 

17 

-1.8171 9355E-0 1 

♦ 

-1.64875002E-01 

- 1 .64733 1 24E-0 1 

* 

-1.97966129E-01 

* 

-1.67412781E-01 

-1.67404319E-01 

18 

-2.04658289E-01 

* 

-1.726493 1  IE-01 

* 

-1.72218483E-0I 

21 

-2.07844699E-01 

15 

-1.75188825E-01 

* 

-1.74695007E-01 

* 

-2.20742079E-01 

♦ 

-1.78192328E-01 

* 

-1.S0676443E-01 

* 

-2.22041721E-01 

* 

-1.81 309666E-0 1 

17 

-1.S2319221E-01 

. 

-2.67948177E-01 

* 

-1.82017191E-01 

* 

- 1 .33 159305E-0 1 

17 

-1.82776874E-01 

-1.S97S0943E-0I 

* 

- 1 .92666 1 39E-0 1 

-1.95481 1 1  IE-01 

♦ 

-1.98326255E-01 

-2.01652740E-01 

* 

-2.02233301E-01 

* 

-2.0224943“E-01 

* 

-2.20500814E-01 

* 

-2.06222481E-C1 

-2.36270762E-01 

* 

-2.340S4842E-0I 

-3. 6246301 8  E-0 1 

* 

-2.6509802SE-01 

* 

-3.88173344E-01 

* 

-3.1  1342750E-O1 

j _ L 
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Table  VI 

Least  stable  eigenvalues  for  a 


Mode  Number 

Eigenvalue  c, 

1 

+  0.00373967 

2 

-  0.03516728 

3 

-  0.03518658 

4 

-  0.05089873 

5 

-  0.06320150 

6 

-  0.06325157 

7 

-  0.09122274 

8 

-  0.09131286 

9 

-  0.11923285 

10 

-  0.11937073 

11 

-  0.12450198 

12 

-  0.13822652 

13 

-  0.1472339 

14 

-  0.1474256 

15 

-  0.1752287 

16 

-  0.1754781 

17 

-  0.1828219 

18 

-  0.203221 

19 

-  0.203529 

20 

-  0.206465 

21 

-  0.208731 

22 

-  0.23119 

23 

-  0.23159 

24 

-  0.23882 

25 

-  0.25872 

26 

-  0.25969 

27 

-  0.25988 

28 

-  0.26511 

29 

-  0.26716 

30 

-  0.28551 

31 

-  0.28663 

32 

-  0.28765 

